source("http://www.statpower.net/Content/312/R Stuff/Steiger R Library Functions.txt")
# 1-way n=3 example  data are 1,2,3   4,5,6   7,8,9
x <- matrix(1:9,9,1)
D1 <- matrix(c(rep(1,3),rep(0,9),rep(1,3),rep(0,9),rep(1,3)),9,3)
H <- t(O(3))
compute.F <- function(x,D1,H){
MS.hypothesis <- t(x) %*% D1 %*% solve(t(D1)%*%D1) %*% H %*% 
solve(t(H) %*% solve(t(D1) %*% D1) %*%H) %*% t(H)%*%
  solve(t(D1)%*%D1) %*% t(D1) %*% x / dim(H)[2]   
MS.error <- t(x) %*% Q(D1) %*% x / (length(x)-dim(H)[1])
F <- MS.hypothesis/MS.error
return(list(F=F[1],df.numerator=dim(H)[2],df.denominator=length(x)-dim(H)[1]))
}

compute.F(x,D1,H)
# 2x2 ANOVA
# 1,2,3  4,5,6
# 7,8,9  13,14,15

#Test row effect
D1 <- matrix(c(rep(1,3),rep(0,12),rep(1,3),rep(0,12),rep(1,3),rep(0,12),rep(1,3)),12,4)
x <- matrix(c(1:9,13:15),12,1)
H <- t(O(2) %x% t(One(2)))
compute.F(x,D1,H)


data <- matrix(c(1,1,1,
                 1,1,2,
                 1,1,3,
                 1,2,4,
                 1,2,5,
                 1,2,6,
                 2,1,7,
                 2,1,8,
                 2,1,9,
                 2,2,13,
                 2,2,14,
                 2,2,15),12,3,byrow=TRUE)
colnames(data) <- c("A","B","Score")
data <- data.frame(data)
fit <- lm(Score ~ factor(A)*factor(B),data=data)
anova(fit)